Deciphering Selectivity Mechanism of BRD9 and TAF1(2) toward Inhibitors Based on Multiple Short Molecular Dynamics Simulations and MM-GBSA Calculations

BRD9 and TAF1(2) have been regarded as significant targets of drug design for clinically treating acute myeloid leukemia, malignancies, and inflammatory diseases. In this study, multiple short molecular dynamics simulations combined with the molecular mechanics generalized Born surface area method were employed to investigate the binding selectivity of three ligands, 67B, 67C, and 69G, to BRD9/TAF1(2) with IC50 values of 230/59 nM, 1400/46 nM, and 160/410 nM, respectively. The computed binding free energies from the MM-GBSA method displayed good correlations with that provided by the experimental data. The results indicate that the enthalpic contributions played a critical factor in the selectivity recognition of inhibitors toward BRD9 and TAF1(2), indicating that 67B and 67C could more favorably bind to TAF1(2) than BRD9, while 69G had better selectivity toward BRD9 over TAF1(2). In addition, the residue-based free energy decomposition approach was adopted to calculate the inhibitor–residue interaction spectrum, and the results determined the gatekeeper (Y106 in BRD9 and Y1589 in TAF1(2)) and lipophilic shelf (G43, F44, and F45 in BRD9 and W1526, P1527, and F1528 in TAF1(2)), which could be identified as hotspots for designing efficient selective inhibitors toward BRD9 and TAF1(2). This work is also expected to provide significant theoretical guidance and insightful molecular mechanisms for the rational designs of efficient selective inhibitors targeting BRD9 and TAF1(2).


Introduction
Acute myeloid leukemia (AML), a morphologically, clinically, and genetically heterogeneous disorder caused by mutations in myeloid differentiation and proliferation, has severely imperiled patients' lives around the world [1][2][3]. Drug design is the creative process of finding specific small molecules that can obstruct or enhance the functions of biological targets based on the action mechanism of the drug and target [4][5][6][7]. Designing small molecules that inhibit the activity of bromodomain-containing proteins (BRDs) is a promising therapeutic strategy to treat many kinds of diseases, including cancer, inflammation, and cardiovascular and autoimmune diseases [8][9][10][11]. Bromodomain-containing protein 9 (BRD9) belonging to the BRD family, a major constituent of the SWI/SNF chromatin remodeling complex named non-canonical BAF, has been identified as a novel therapeutic target in AML [12,13]. The sensitivity of AML cells is correlated with the level of BRD9, and AML cells endure terminal differentiation and cycle arrest with the degradation of BRD9 [14]. Meanwhile, the second bromodomain of the human transcription initiation factor TFIID subunit 1 (TAF1 (2)) is overexpressed in a variety of cancers and plays a significant role in AML1-ETO fusion gene expression [15]. Furthermore, multiple reports indicate the key roles of TAF1 (2) in AML and provide a new theoretical structural framework to develop direct-acting small molecule inhibitors of TAF1(2) as prospective inflammation pathophysiology and cancer therapeutics [16][17][18][19][20]. framework to develop direct-acting small molecule inhibitors of TAF1(2) as pro inflammation pathophysiology and cancer therapeutics [16][17][18][19][20].
The BRD9 and TAF1 (2), as shown in Figure 1, structurally share four left-ha helices (αA, αB, αC, αZ) constituting up-and-down four-helix bundles, which f loops between helices αA and αZ (ZA loop) and αB and αC (BC loop), respectiv These two loops feature a hydrophobic pocket, which frequently generates two s drogen bonds, one is between the amide and asparagine at the top of the BC BRD9/TAF1 (2), and the other is water-mediated and formed with the tyrosine sit the ZA loop of BRD9/TAF1 (2) [22]. In addition, the hydrophobic binding p BRD9/TAF1 (2) consists of the gatekeeper at the head of the αC helix and the li shelf with the first three amino acid residues of the ZA loop [23]. In BRD9, three r G43, F44, and F45, adjacent to the ZA channel constitute the lipophilic shelf, and the beginning of the αC helix is the gatekeeper. Meanwhile, the lipophilic shelf of is composed of three residues, W1526, P1527, and F1528, and residue Y1589 is r as the gatekeeper. Although the tyrosine and asparagine residues are heavily co in the vast majority of bromodomain proteins, there are apparent conformational in the gatekeeper, lipophilic shelf, and ZA channel residues. Based on the importa roles of BRD9 and TAF1 (2) in drug design toward human cancers, Crawford and ers solved the crystal structures of 67B-and 69G-bound BRD9, as well as 67B-a associated TAF1(2) [22]. Despite the highly similar structures shared by 67B, 67C, ( Figure 2C-E), three inhibitors have different binding affinities to BRD9 and TAF1 IC50 values of 230/59 nM, 1400/46 nM, and 160/410 nM for BRD9/TAF1(2), resp Therefore, it is essential to clarify the binding selectivity of inhibitors to BRD9 and and the conformational changes of the two proteins caused by inhibitor binding signing drugs for anti-cancer treatment. Tertiary structure superimposition of the BRD9 (PDB code 5I7X) and TAF1(2) (P 5I29) bromodomains, in which BRD9 is displayed in green and TAF1 (2) in magenta wit modes. Notable residues include the gatekeeper (Y106 in BRD9 and Y1589 in TAF1(2)), the shelf adjacent to the ZA channel (G43, F44, and F45 in BRD9), as well as W1526, P1527, a in TAF1(2). Tertiary structure superimposition of the BRD9 (PDB code 5I7X) and TAF1(2) (PDB code 5I29) bromodomains, in which BRD9 is displayed in green and TAF1 (2) in magenta with cartoon modes. Notable residues include the gatekeeper (Y106 in BRD9 and Y1589 in TAF1(2)), the lipophilic shelf adjacent to the ZA channel (G43, F44, and F45 in BRD9), as well as W1526, P1527, and F1528 in TAF1 (2).  (2), and 67B in complex with BRD9 (PDB code 5I7X) and 67B bound to TAF1(2) (PDB 5I29), in which BRD9 is displayed in green and TAF1 (2) in magenta with cartoon modes, (B) Bin pockets of BRD9 and TAF1 (2) are shown in the surface modes and inhibitor in the stick mode E) correspond to the structures of 67B, 67C, and 69G, respectively, from which inhibitors ar picted in line forms.
Until now, various simulation approaches, including conventional molecular namics (cMD) [24][25][26][27][28][29], multiple short molecular dynamics (MSMD) [30], accelerated lecular dynamics (aMD) [31][32][33], Gaussian accelerated molecular dynamics (GaMD) 38], and multiple replica Gaussian accelerated molecular dynamics (MR-GaMD) sim tions [39][40][41] have been used to investigate the conformational alterations and bin mechanisms of receptors due to ligand associations and residue mutations [42][43][44]. eral intensive molecular dynamics studies have successfully deciphered the molec mechanism regarding the binding selectivity of inhibitors toward homological pro with highly similar tertiary structures [45][46][47][48][49][50][51][52]. In fact, numerous molecular dynamics ulation works were also conducted to research the binding modes of different ligan BRD9 and TAF1 (2). For instance, Wang et al. combined molecular dynamics and bin free energy predictions to determine the binding mechanism of three small molecule ands, 5SW, 5U2, and 5U6, toward BRD9 and identify the hot interaction spots of the tein with inhibitors [53]. Liu's group applied MSMD simulations and binding affinity culations to the investigation of the binding difference of inhibitors to different BRD ilies, and their results provided useful information for clarifying the selective mechan of BRD families [54][55][56]. Song et al. estimated the binding free energies and energetic tributions of individual residues of four pyridinone-like scaffold inhibitors compl with BRD9 based on cMD simulations, and their results suggested that the aromatic and dimethoxyphenyl structure combined into a pyridinone scaffold effectively enha the BRD9 binding affinity [57]. Magno et al. performed twenty-four independent MD ulations and free energy profile analyses to investigate the spontaneous and rever  (2), and 67B in complex with BRD9 (PDB code 5I7X) and 67B bound to TAF1(2) (PDB code 5I29), in which BRD9 is displayed in green and TAF1(2) in magenta with cartoon modes, (B) Binding pockets of BRD9 and TAF1(2) are shown in the surface modes and inhibitor in the stick mode, (C-E) correspond to the structures of 67B, 67C, and 69G, respectively, from which inhibitors are depicted in line forms.
Until now, various simulation approaches, including conventional molecular dynamics (cMD) [24][25][26][27][28][29], multiple short molecular dynamics (MSMD) [30], accelerated molecular dynamics (aMD) [31][32][33], Gaussian accelerated molecular dynamics (GaMD) [34][35][36][37][38], and multiple replica Gaussian accelerated molecular dynamics (MR-GaMD) simulations [39][40][41] have been used to investigate the conformational alterations and binding mechanisms of receptors due to ligand associations and residue mutations [42][43][44]. Several intensive molecular dynamics studies have successfully deciphered the molecular mechanism regarding the binding selectivity of inhibitors toward homological proteins with highly similar tertiary structures [45][46][47][48][49][50][51][52]. In fact, numerous molecular dynamics simulation works were also conducted to research the binding modes of different ligands to BRD9 and TAF1 (2). For instance, Wang et al. combined molecular dynamics and binding free energy predictions to determine the binding mechanism of three small molecule ligands, 5SW, 5U2, and 5U6, toward BRD9 and identify the hot interaction spots of the protein with inhibitors [53]. Liu's group applied MSMD simulations and binding affinity calculations to the investigation of the binding difference of inhibitors to different BRD families, and their results provided useful information for clarifying the selective mechanisms of BRD families [54][55][56]. Song et al. estimated the binding free energies and energetic contributions of individual residues of four pyridinone-like scaffold inhibitors complexed with BRD9 based on cMD simulations, and their results suggested that the aromatic ring and dimethoxyphenyl structure combined into a pyridinone scaffold effectively enhanced the BRD9 binding affinity [57]. Magno et al. performed twenty-four independent MD simulations and free energy profile analyses to investigate the spontaneous and reversible binding of acetylated lysine to TAF1 (2), and the obtained dynamical information indicated that hydrogen bond interactions stabilized the two main binding modes of TAF1(2) [58]. Thus, it is important to decipher the binding mode and selectivity mechanism of inhibitors at the atomic levels for the development of potent small-molecule inhibitors targeting BRD9 and TAF1 (2).
In the current study, in order to illustrate the selective mechanism of BRD9 and TAF1(2), three small-molecule inhibitors, 67B, 67C, and 69G, were chosen to determine their binding selectivity for BRD9 and TAF1(2) by performing MSMD simulations and binding free energy computations. Encouragingly, MSMD simulations can extract more rational conformational sampling than cMD simulations, which has been verified in previous works [59][60][61]. Furthermore, principal component analysis (PCA) [62,63], dynamics crosscorrelation maps (DCCMs), and free energy landscapes (FELs) were combined to investigate the conformational variations and internal dynamics of BRD9 and TAF1(2) caused by inhibitor associations. To contrast the structural differences between BRD9 and TAF1(2), the PyMOL software was utilized to align the two complexes, and their structures are displayed in Figure 2A. As shown in Figure 2B, the binding pockets of BRD9 and TAF1 (2) were drawn in surface forms, while the ligand was depicted in stick form. The structures of three small-molecule inhibitors. 67B, 67C, and 69G, are shown in line forms in Figure 2C-E, respectively. The three inhibitors 67B, 67C, and 69G share similar structures, except for the red rectangle. It was conducive to determine the impact of the structural variations in 67B, 67C, and 69G on the binding selectivity of BRD9 and TAF1(2) for the design of small-molecule inhibitors associated with bromodomain proteins. In this work, MSMD simulations and several analysis methods, such as PCA, inhibitor-residue interactions, calculations of DCCMs, analysis of FELs, hydrogen bonding interactions (HBIs), and hydrophobic interactions (HIs), were integrated to carry out this study. In the meantime, this study is also expected to provide rational information for the design of small-molecule inhibitors toward BRD9 and TAF1(2).

Structural Flexibilities and Fluctuations of BRD9 and TAF1(2)
In order to execute reliable and rational conformational samplings of BRD9 and TAF1(2), a total of 1.2 µs MSMD simulations, including three individual cMD simulations of 400 ns, were conducted on the apo BRD9 and TAF1(2), as well as the six inhibitor-BRD9 and inhibitor-TAF1(2) systems with three inhibitors, 67B, 67C, and 69G. To assess the overall stability of the MSMD simulations, the fluctuations in the root-mean-squaredeviations (RMSDs) of backbone atoms in BRD9 and TAF1(2) relative to the corresponding initial conformations of the six complexes over time were calculated and are depicted in Figure S1 see in Supplementary Materials. On the whole, the structural fluctuation range of the apo BRD9 was higher than that of the inhibitor-BED9 complexes while the structural fluctuation extent of the apo TAF1(2) was similar to that of the inhibitor-TAF1(2) complexes ( Figure S1A,B). The structural variations of three replicas of 67B-, 67C-, and 69Gassociated BRD9/TAF1(2) were convergent after 100 ns of cMD simulations ( Figure S1C-H). Hence, the stable portions (100-400 ns) from three independent cMD simulations were concatenated together to create a single integrated trajectory (SIT) of 900 ns for each complex, which was used to execute all computations and dynamics analyses.
To further determine the structural flexibilities of the BRD9 and TAF1(2) induced by inhibitor associations, the root-mean-square fluctuations (RMSFs) of the C α atoms in BRD9 and TAF1(2) were calculated based on the SIT (Figure 3). According to the comparison, BRD9 and TAF1(2) produced similar RMSFs fluctuations, suggesting that BRD9 and TAF1(2) embodied common rigid and flexible domains. For the BRD9-related systems, the binding of three inhibitors weakened the structural flexibility of BRD9, especially for the ZA-loop ( Figure 3A,B). However, for the TAF1(2)-related systems, the binding of 67B and 67C slightly reduced the structural flexibility of TAF1(2), while the presence of 69G strengthened that of TAF1(2), particularly for the ZA-loop ( Figure 3C,D). Obvious structural alterations mainly existed in four regions, including L1 (residues 36-45 for BRD9 and 1519-1528 for TAF1(2)), L2 (residues 49-67 for BRD9 and 1532-1550 for TAF1(2)), L3 (residues 74-85 for BRD9 and 1557-1568 for TAF1(2)), and L4 (residues 99-110 for BRD9 and 1582-1593 for TAF1 (2)). The lipophilic shelf contained three residues (G43, F44, and F45 in BRD9, and W1526, P1527, and F1528 in TAF1(2)) at the end of the region L1, while the gatekeeper included one residue (Y106 in BRD9, and Y1589 in TAF1(2)) at region L4. For BRD9, the inhibitor associations induced evident alterations in the structural flexibility at regions L1 and L2 ( Figure 3A); however, the bindings of inhibitors with TAF1(2) only yielded an obvious influence on region L2 ( Figure 3C). These alterations in RMSFs indicated that the structural flexibility of BRD9 was higher than that of TAF1 (2). The RMSF values of the 67B-, 67C-, and 69G-BRD9 compounds in region L3 and the corresponding ones in region L4 of TAF1(2) were lower than 1.0 Å, demonstrating that the two regions were rigid. However, the RMSFs values of region L2 in BRD9 and TAF1(2) were above 1.5 Å, suggesting that region L2 was flexible. Owing to 67B and 67C binding, the main portions of RMSFs in BRD9 were obviously larger than the corresponding ones in TAF1 (2). In the bound state of 69G, the flexibility of three regions, L1, L2, and L3, in BRD9 was lower than those of TAF1(2), while the flexibility of region L4 in BRD9 was slightly higher than that of TAF1 (2). The results signify that several residues of the four aforementioned regions were potential central factors driving the selective binding of inhibitors toward these two proteins.
In the bound state of 69G, the flexibility of three regions, L1, L2, and L3, in BRD9 was lower than those of TAF1(2), while the flexibility of region L4 in BRD9 was slightly higher than that of TAF1 (2). The results signify that several residues of the four aforementioned regions were potential central factors driving the selective binding of inhibitors toward these two proteins.  As shown in Figure 3B, G43 was the first residue of ZA channel, Y106 was the gatekeeper of BRD9, W1526 was the first residue of ZA channel, and Y1589 was the corresponding gatekeeper in TAF1 (2). The distances between the first residue of the ZA channel and the gatekeeper (Y106-G43 in BRD9 or Y1589-W1526 in TAF1(2)) were determined from the SIT, and the corresponding frequency distributions are displayed in Figure 4A,C,E.
Moreover, the frequency distributions of the Chi dihedral angle of the side chain (Y106 in BRD9 and Y1589 in TAF1 (2)) are depicted in Figure 4B,D,F. As shown in Figure 4A, the higher peak values of the Y106 C α -G43 C α distance of the BRD9-67B complex were distributed range between 11.1 and 12.0 Å, while the higher peak of the Y1589 Cα-W1526 Cα distance of the TAF1(2)-67B complex was at 10.2 Å. The above results demonstrate that the distance between the gatekeeper and the ZA channel in BRD9 was slightly larger than that in TAF1 (2), which suggested that the hot spot site volume of TAF1(2) was smaller than that of BRD9. The peak values of the Chi dihedral angle in the side chain of Y106 for the 67B-, 67C-, and 69G-BRD9 complexes were 313.3 • , 310.0 • , and 308.0 • , respectively, while the ones for the 67B-, 67C-, and 69G-TAF1(2) complexes were 297.8 • , 299.7 • , and 301.8 • , separately. The results show that the hot spots of BRD9 had higher flexibility than those of TAF1(2), which is consistent with the above fluctuations in the RMSFs.
As shown in Figure 3B, G43 was the first residue of ZA channel, Y106 was the gatekeeper of BRD9, W1526 was the first residue of ZA channel, and Y1589 was the corresponding gatekeeper in TAF1 (2). The distances between the first residue of the ZA channel and the gatekeeper (Y106-G43 in BRD9 or Y1589-W1526 in TAF1(2)) were determined from the SIT, and the corresponding frequency distributions are displayed in Figure  4A,C,E. Moreover, the frequency distributions of the Chi dihedral angle of the side chain (Y106 in BRD9 and Y1589 in TAF1 (2)) are depicted in Figures 4B, D, and F. As shown in Figure 4A, the higher peak values of the Y106 Cα-G43 Cα distance of the BRD9-67B complex were distributed range between 11.1 and 12.0 Å, while the higher peak of the Y1589 Cα-W1526 Cα distance of the TAF1(2)-67B complex was at 10.2 Å . The above results demonstrate that the distance between the gatekeeper and the ZA channel in BRD9 was slightly larger than that in TAF1 (2), which suggested that the hot spot site volume of TAF1(2) was smaller than that of BRD9. The peak values of the Chi dihedral angle in the side chain of Y106 for the 67B-, 67C-, and 69G-BRD9 complexes were 313.3°, 310.0°, and 308.0°, respectively, while the ones for the 67B-, 67C-, and 69G-TAF1(2) complexes were 297.8°, 299.7°, and 301.8°, separately. The results show that the hot spots of BRD9 had higher flexibility than those of TAF1(2), which is consistent with the above fluctuations in the RMSFs.

Internal Dynamics of BRD9 and TAF1(2)
To explore the changes in the internal dynamics of BRD9 and TAF1(2) due to inhibitor associations, the cross-correlation coefficients were calculated by using the C α atomic coordinates recorded in the SIT, and cross-correlation maps are displayed in Figure 5. According to the color-decoded patterns, the negative regions (dark blue and plain blue) indicated extremely anticorrelated (AC) movements, while the highly positive regions (red and yellow) were related to strongly positively correlated (PC) motions between specific residues. As shown in Figure 5, the binding of 67B, 67C, and 69G produced evident influences on the motion modes of BRD9 and TAF1 (2).

Internal Dynamics of BRD9 and TAF1(2)
To explore the changes in the internal dynamics of BRD9 and TA tor associations, the cross-correlation coefficients were calculated by coordinates recorded in the SIT, and cross-correlation maps are displa cording to the color-decoded patterns, the negative regions (dark b indicated extremely anticorrelated (AC) movements, while the high (red and yellow) were related to strongly positively correlated (PC) m cific residues. As shown in Figure 5, the binding of 67B, 67C, and 69 influences on the motion modes of BRD9 and TAF1(2). For BRD9 associated with 67B, 67C, and 69G ( Figure 5A,C,E), reg generated significant AC motions. Compared with BRD9 complexed 69G, the binding of 67B, 67C, and 69G to TAF1(2) reduced the AC mo R2, and R3 ( Figure 5B,D,F). The results demonstrate that the binding tors yielded distinct influences on the internal dynamics of BRD9 and For BRD9 associated with 67B, 67C, and 69G ( Figure 5A,C,E), regions R1, R2, and R3 generated significant AC motions. Compared with BRD9 complexed with 67B, 67C, and 69G, the binding of 67B, 67C, and 69G to TAF1(2) reduced the AC motions in regions R1, R2, and R3 ( Figure 5B,D,F). The results demonstrate that the bindings of identical inhibitors yielded distinct influences on the internal dynamics of BRD9 and TAF1 (2), which signified that important residues located in R1-R3 of BRD9 and TAF1(2) may have yielded significant hydrophobic and hydrogen bonding interactions with inhibitors and played critical roles in the binding selectivity of ligands toward BRD9 and TAF1 (2).
In addition, principal component analysis (PCA) was conducted to decode the conformational alterations of BRD9 and TAF1(2) due to the associations with 67B, 67C, and 69G, respectively, and the function of forty eigenvalues stemming from the diagonalization of covariance matrix versus the related eigenvector indexes is displayed in Figure 6. The first few larger eigenvalues represent the primarily collective motions of the structural domain in these two proteins. The first six eigenvalues accounted for 79.07%, 72.72%, and 60.36% of the total movements for the 67B-, 67C-, and 69G-BRD9 complexes and 70.37%, 75.08%, and 75.07% of the total motions of the 67B-, 67C-, and 69G-TAF1(2) complexes, respectively.
Molecules 2023, 28, x FOR PEER REVIEW  8 of 27 69G, respectively, and the function of forty eigenvalues stemming from the diagonalization of covariance matrix versus the related eigenvector indexes is displayed in Figure 6. The first few larger eigenvalues represent the primarily collective motions of the structural domain in these two proteins. The first six eigenvalues accounted for 79.07%, 72.72%, and 60.36% of the total movements for the 67B-, 67C-, and 69G-BRD9 complexes and 70.37%, 75.08%, and 75.07% of the total motions of the 67B-, 67C-, and 69G-TAF1(2) complexes, respectively. To gain more insight into the alterations in motional modes between BRD9 and TAF1(2) induced by inhibitor associations, the first eigenvectors of six systems were visualized in six porcupine plots ( Figure S3). The direction of the arrow reflects the direction of the movements and the length of the arrow represents the strength of the motions. In contrast with 67B-, 67C-, and 69G-bound BRD9, the bindings of these inhibitors not only altered the movement directions of the L2, L3, and L4 in TAF1 (2), but also altered the movement amplitude of these three loops. Furthermore, the αZ helix of the 67B-BRD9 complex moved toward the left and down ( Figure S3A), while that of the 67B-TAF1 (2) was altered toward the right and up ( Figure S3B). The αZ helix of the 67C-BRD9 complex moved toward the right ( Figure S3C), while that of the 67C-TAF1(2) was transformed toward the left and down ( Figure S3D). The αZ helix of the 69G-BRD9 complex moved upward ( Figure S3E), but that of the 69G-TAF1(2) was altered toward the left and down ( Figure S3F). The above discussions indicate that the conformational alterations of BRD9 and TAF1(2) extracted from MSMD simulations probably led to the distinct binding selectivities of 67B, 67C, and 69G toward BRD9 and TAF1(2). To gain more insight into the alterations in motional modes between BRD9 and TAF1(2) induced by inhibitor associations, the first eigenvectors of six systems were visualized in six porcupine plots ( Figure S3). The direction of the arrow reflects the direction of the movements and the length of the arrow represents the strength of the motions. In contrast with 67B-, 67C-, and 69G-bound BRD9, the bindings of these inhibitors not only altered the movement directions of the L2, L3, and L4 in TAF1 (2), but also altered the movement amplitude of these three loops. Furthermore, the αZ helix of the 67B-BRD9 complex moved toward the left and down ( Figure S3A), while that of the 67B-TAF1(2) was altered toward the right and up ( Figure S3B). The αZ helix of the 67C-BRD9 complex moved toward the right ( Figure S3C), while that of the 67C-TAF1(2) was transformed toward the left and down ( Figure S3D). The αZ helix of the 69G-BRD9 complex moved upward ( Figure S3E), but that of the 69G-TAF1(2) was altered toward the left and down ( Figure S3F). The above discussions indicate that the conformational alterations of BRD9 and TAF1(2) extracted from MSMD simulations probably led to the distinct binding selectivities of 67B, 67C, and 69G toward BRD9 and TAF1(2).

Binding Ability of Inhibitors to BRD9 and TAF1(2)
To further evaluate the variance in the binding abilities of 67B, 67C, and 69G to BRD9 and TAF1(2), the MM-GBSA approach was employed to calculate the binding free energies (BFEs) of three ligands to BRD9 and TAF1(2) by using 300 conformational structures withdrawn from the 900 ns SIT with a time step of 3 ns. Fifty structural frames were taken from the above 300 conformations at an interval of 6 conformations to calculate the contributions of entropy (−T∆S) to the binding associations through the normal mode analysis approach. All energetic components resulting from the MM-GBSA calculations are listed in Table 1. The ranks of the BFEs of 67B-BRD9/TAF1(2), 67C-BRD9/TAF1(2), and 69G-BRD9/TAF1(2) were consistent with those of the experimental values, which demonstrated that the computed BFEs were reliable and rational. The energies with positive values provided unfavorable factors for inhibitor associations, while the negative components contributed favorable forces for inhibitor bindings. As shown in Table 1, for BRD9 and TAF1(2) bound by 67B, 67C, and 69G, the negative electrostatic interaction energies (∆E ele ) were absolutely overwhelmed by positive polar solvation energies (∆G gb ) to form unfavorable terms (∆G ele+gb ) for inhibitor associations. The unfavorable factors of the entropy changes (−T∆S) also weakened the binding associations of 67B, 67C, and 69G to BRD9 and TAF1 (2). Meanwhile, the negative values of the nonpolar solvation energies (∆G nonpol ) and Van der Waals interactions (∆E vdW ) produced favorable factors (∆E vdW+nonpol ) for the inhibitor-BRD9/TAF1(2) bindings. As shown in Table 1, the value of ∆G ele+gb for 67B-TAF1(2) was reduced by 0.90 kcal/mol relative to that of 67B-TAF1(2), and the favorable term ∆E vdW+nonpol of 67B-TAF1(2) was enhanced by 2.77 kcal/mol in comparison with that of 67B-BRD9, which led to an increase of 3.67 kcal/mol in the enthalpy changes of the 67B-TAF1(2) compared with that of 67B-BRD9. Moreover, the value of −T∆S for 67B-TAF1(2) was strengthened by 0.94 kcal/mol relative to that of 67B-BRD9. Overall, the binding affinity of 67B-TAF1(2) was enhanced by 2.73 kcal/mol, suggesting that 67B generated a stronger association with TAF1(2) than with BRD9. In the case of 67C, an increase of 3.41 kcal/mol in the ∆E vdW+nonpol of the 67C-TAF1(2) complex compared with that of 67C-BRD9 led to an increase in the enthalpy changes of the 67C-TAF1(2) complex compared with that of 67C-BRD9. The entropy changes −T∆S of the 67C-TAF1(2) complex decreased by 1.94 kcal/mol relative to that of 67C-BRD9, which finally resulted in an increase of 4.99 kcal/mol in the binding free energy of the 67C-TAF1(2) complex relative to that of 67C-BRD9. Therefore, the binding affinity of 67C to TAF1(2) was higher than that for BRD9. For ligand 69G, the value of ∆G ele+gb of the 69G-TAF1(2) complex was decreased by 0.6 kcal/mol compared with that of the 69G-BRD9 complex, and the negative component ∆E vdW+nonpol of the 69G-TAF1(2) complex was decreased by 2.94 kcal/mol relative that of 69G-BRD9, which resulted in an overall decrease of 2.34 kcal/mol in the binding enthalpy of the 69G-TAF1(2) complex relative to that of 69G-BRD9. Furthermore, the −T∆S of the 69G-TAF1(2) complex was decreased by 0.57 kcal/mol compared with that of 69G-BRD9. In view of the above two factors, the binding ability of 69G to BRD9 was increased by 1.77 kcal/mol relative to TAF1(2), suggesting that 69G produces more favorable binding to BRD9 than TAF1(2).

Alterations in the Free Energy Landscapes of BRD9 and TAF1(2) Caused by Inhibitor Bindings
The Free energy landscapes (FELs) can effectively represent various free energy states relative to the conformational alterations of proteins due to changes in the binding environment [64][65][66]. To study the influences of small molecular inhibitors' associations on the conformational alterations of BRD9 and TAF1(2), projections of the SIT onto the first two principal components arising from the diagonalization of the covariance matrix were used as reaction coordinates to construct the FELs of BRD9 and TAF1(2). The results and structures are displayed in Figure 12, Figure 13 and Figure S6, in which the symbols V1, V2, V3, and V4, and the red points denote different energy valleys recognized by MSMD simulations. As shown in these figures, the associations of three inhibitors, 67B, 67C, and 69G, with BRD9 and TAF1(2) evidently affected the FELs and generated conformational alterations. The MSMD simulations captured two distinct energy valleys, V1 and V2, and according to the color bar, the typical structure V2 was located at a deeper potential basin than structure V1 (Figure 12A), indicating that the conformations of the 67B-BRD9 complex were primarily distributed at two energetic spaces. The structures of 67B in two typical conformations V1 and V2 of the 67B-associated BRD9 are aligned together in Figure 12B. The superimpositions of two representative structures corresponding to energy valleys V1 and V2 reveal that three domains, L1, L2, and L4, significantly deviated from each other ( Figure 12B), suggesting that 67B generates an evident effect on the structural flexibility of the 67B-BRD9 complex. The results showed that 67B had two various binding poses in BRD9, which significantly influenced the association of 67B to BRD9. Regarding the 67B-associated TAF1(2), two distinct energy valleys, V1 and V2, were identified by using entire MSMD simulations, and the color bar indicated that the depth of energy valley V1 was deeper than that of energy valley V2 in Figure 13A, indicating that the structures of the 67B-TAF1(2) complex were mainly distributed at two energetic spaces. The alignments of two conformations located at energy valleys V1 and V2 demonstrated that domains L1 and L2 produced evident deviations from each other ( Figure 13B), which possibly generated significant effects on the association of 67B with TAF1(2). The conformational superimpositions of two representative structures of the 67B-TAF1(2) complex in energy valleys V1 and V2 denoted that the association poses of 67B in TAF1(2) yielded primarily parallel slides from each other, which exerted a specific influence on the binding of 67B to TAF1(2).   The MSMD simulations detected four primary energy valleys, V1, V2, V3, and V4, in the 67C-associated BRD9; in accordance with the color bar, the depth of energy valley V1 was the deepest in these four energy valleys ( Figure 12C), suggesting that the structures of the 67C-BRD9 complex were chiefly populated at four energetic spaces. Four typical structures located at valleys V1, V2, V3, and V4 are superposed in Figure 12D, and the results demonstrate that loops L1 and L2 significantly diverged from each other, indicating that these two loops displayed a considerable structural flexibility and played an essential role in the association of 67C with BRD9. As shown in the conformational alignment of 67C in the typical structures V1, V2, V3, and V4 ( Figure 12D), 67C had four distinct association poses and produced large deviations. For the 67C-TAF1(2) complex, three energy valleys, V1, V2, and V3, were recognized by the entire MSMD simulations ( Figure  13C), demonstrating that the 67C-TAF1(2) complex was distributed at three conformational subspaces. The alignments of three typical conformations situated at energy basins V1, V2, and V3 denote that domains L1 and L2 yielded evident distortions from each other The MSMD simulations detected four primary energy valleys, V1, V2, V3, and V4, in the 67C-associated BRD9; in accordance with the color bar, the depth of energy valley V1 was the deepest in these four energy valleys ( Figure 12C), suggesting that the structures of the 67C-BRD9 complex were chiefly populated at four energetic spaces. Four typical structures located at valleys V1, V2, V3, and V4 are superposed in Figure 12D, and the results demonstrate that loops L1 and L2 significantly diverged from each other, indicating that these two loops displayed a considerable structural flexibility and played an essential role in the association of 67C with BRD9. As shown in the conformational alignment of 67C in the typical structures V1, V2, V3, and V4 ( Figure 12D), 67C had four distinct association poses and produced large deviations. For the 67C-TAF1(2) complex, three energy valleys, V1, V2, and V3, were recognized by the entire MSMD simulations ( Figure 13C), demonstrating that the 67C-TAF1(2) complex was distributed at three conformational subspaces. The alignments of three typical conformations situated at energy basins V1, V2, and V3 denote that domains L1 and L2 yielded evident distortions from each other ( Figure 13D).
Active regions L1 and L2 yielded obvious distortions among three typical conformations, which probably brought significant effects on the poses of 67C. This result is supported by the apparent structural slides and torsions of 67C shown in Figure 13D. Therefore, the conformational alterations of L1 and L2, and slides and torsions of 67C certainly exerted evident influences on the binding selectivity of 67C to BRD9 and TAF1(2).
2.5.3. The 69G-Bound BRD9 over the 69G-Associated TAF1 (2) Two energy valleys, V1 and V2, were captured by the MSMD simulations of the 69G-BRD9 complex and on the basis of the color bar; structure V2 was located at a deeper valley than structure V1 ( Figure 12E), indicating that the 69G-BRD9 complex occupied two primary energetic spaces. The superimpositions of two representative conformations located at energy valleys V1 and V2 suggested that loop L1 had higher structural flexibility and apparently diverged from each other in the 69G-BRD9 system ( Figure 12F). The conformations of 69G in two representative basins, V1 and V2, were superimposed ( Figure 12F) and the results indicate that 69G yielded apparent sliding, which indicates a significant influence on the associations of 69G with BRD9. For the 69G-bound TAF1(2), four energy valleys, V1, V2, V3, and V4, were identified by the whole MSMD simulation, and in accordance with the color bar, the depth of the energy valleys in state V3 was deeper than that of the other three energy valleys, V1, V2, and V3 ( Figure 13E), indicating that the conformations of the 69G-TAF1(2) were primarily clustered into four energetic spaces. The alignments of the 69G-TAF1(2) complex located at energy valleys V1, V2, V3, and V4 indicated that loops L1 and L2 produced obvious deviations from each other ( Figure 13F). The conformations of 69G in four typical structures of the 69G-TFA1(2) complex were superimposed ( Figure 13F), and the results demonstrated that 69G contained four distinct association poses and produced evident distortions, which obviously affected the association of 69G with TAF1(2). Therefore, the changes in the orientations and conformations explicitly influenced the binding selectivity of ligand 69G toward BRD9 and TAF1(2).

Modeling Simulated Systems
The initial configurations of the crystal structures were obtained from the Protein Data Bank (PDB): 5I7X and 5I7Y corresponded to the 67B-and 69G-BRD9 complexes, respectively, while 5I29 and 5I1Q corresponded to the 67B-and 67C-TAF1(2) complexes, respectively [22]. Meanwhile, as the crystal structures of the 67C-BRD9 and 69G-TAF1(2) complexes were unavailable, the crystal structures 517Y were superimposed with 5I1Q to generate the configuration of the 670C-BRD9 compound by deleting 69G and TAF1(2) through the PyMol software (https://www.pymol.org accessed on 5 February 2023). Similarly, the crystal structures of the 69G-BRD9 (5I7Y) and 67C-TAF1(2) (5I1Q) complexes were aligned together to produce the initial atomic coordinates of the 69G-TAF1(2) compound by eliminating BRD9 and 67C. Owing to the difference in the length of the residues from BRD9 and TAF1(2), residues 22-102 in BRD9 and residues 1504-1604 in TAF1(2) were used as the starting crystal structures of the MSMD simulations. All missing hydrogen atoms in the crystal structures were added to their corresponding heavy atoms with the Leap module of Amber 20.0 [67,68]. The ff19SB force field [69] and TIP3P model [70] were used to generate the force field parameters of proteins BRD9 and TAF1(2), water molecules, and counter ions. The configurations of three inhibitors, 67B, 67C, and 69G, were optimized using the semi-empirical AM1 approach, and then the atomic BCC charges of 67B, 67C, and 69G were produced by using the Antechamber module in Amber 20.0. The general AMBER force field (GAFF) was applied to generate the force field parameters of three inhibitors, 67B, 67C, and 69G [71]. Four chlorine ions (Cl − ) and eleven sodium ions (Na + ) were deposited around the ligand-associated BRD9 and TAF1(2) to generate six neutral simulated systems in the salt environment of 0.15 M NaCl. In addition, octahedral periodic water boxes with 12.0 Å using the TIP3P model were utilized to solvate the ligand-BRD9 or TAF1(2) complexes, and the number of water molecules was about 7600.

Multiple Short Molecular Dynamics (MSMD) Simulations
pmemd.cuda embedded in Amber 20 was used to implement all of the multiple shortmolecular dynamics simulations [72]. To relax each system, a 2500-step steepest descent minimization, another 2500-step conjugate gradient minimization, a 2 ns soft heating process of 0 to 300 K under constant number, volume, and temperature (NVT) condition, and then a 2 ns equilibrium process of 310 K under the constant number, pressure, and temperature (NPT) condition were further carried out. Ultimately, three independent 400 ns cMD simulations were executed with periodic boundary conditions and the particle mesh Ewald (PME) method at a constant temperature (300K) and pressure (1 bar) to relax each simulated complex [73,74]. During the whole MSDM simulations in the present work, the chemical bonds between hydrogen and heavy atoms were constrained with the SHAKE numerical integration algorithm [75]. The temperatures of the simulated complexes were controlled by utilizing the Langevin equation with a collision frequency of 2.0 ps −1 with a mollified impulse approach for the Newtonian molecular dynamics [76]. A cutoff value of 10 Å was applied to conduct the estimations of the electrostatic interactions with the PME approach and computations of the van der Waals interactions. In the present work, 1.2 µs MSMD simulations, including three individual cMD simulations of 400 ns, were conducted for the ligand-BRD9/TAF1(2) complexes. To execute the post-process investigation, the equilibrium parts from three independent MSMD trajectories were linked into an SIT. The PCA and DCCMs were executed based on the SIT with the CPPTRAJ module in Amber [77], and the details were introduced in our previous works [78,79].

MM-GBSA Free Energy Computations and Decomposition
Among the methods of binding affinity prediction, although the thermodynamics integration [80][81][82] and free energy perturbation [83][84][85][86][87][88][89] methods can provide more accurate results, these two methods are very expensive in terms of computation time. It is highly important to achieve a good trade-off between accuracy and efficiency in calculations of inhibitor-target binding free energies for drug development, which has been discussed by Rizzuti et al. in their work [90]. Recently, empirical equation-based MM-GBSA and molecular mechanics Poisson Boltzman surface area (MM-PBSA) methods have extensively been used to estimate the binding free energy for numerous ligand-protein and proteinprotein interactions [91][92][93][94][95][96][97]. Furthermore, Hou's group assessed the performance of the MM-GBSA and MM-PBSA approaches by calculating the binding free energies of various biological systems, and their results implied that the MM-GBSA approach could provide more reasonable conclusions [98,99]. Therefore, the MM-GBSA approach was employed to compute the binding free energies of 67B, 67C, and 69G to BRD9 and TAF1(2) with the following Equation (1): where the first two components ∆E ele and ∆E vdW denote the electrostatic and van der Waals interactions of 67B, 67C, and 69G with BRD9 and TAF1(2), respectively, which were estimated by using the FF19SB force field. The terms ∆G gb and ∆G nonpol indicated the polar and nonpolar solvation free energies, respectively. The third element was solved based on the GB-OBCI model developed by Onufriev et al., and the fourth term was computed using the following empirical formula [100]: in which parameters γ and ∆SASA denote the surface tension and the changes in the solvent-accessible surface areas due to inhibitor associations, respectively. Parameters γ and β were to as 0.0072 kcal·mol·Å −2 and 0 kcal·mol −1 in this work, respectively [101].

Conclusions
Clarifying the binding selectivity of inhibitors to BRD9 and TAF1(2) plays an important role in drug targets for AML, human malignancies, and inflammatory disease therapy. This study aimed to decipher the molecular mechanism of the binding selectivity of three ligands, 67B, 67C, and 69G, toward BRD9 and TAF1 (2). MSMD simulations of 1.2 µs, including three independent cMD simulations of 400 ns, were performed on six inhibitor-associated BRD9 and TAF1(2) complexes to decode the binding selectivity of small molecular inhibitors to BRD9 and TAF1 (2). The internal dynamics of BRD9 and TAF1(2) were investigated by means of DCCMs, PCA, and FELs, and the results demonstrated that the associations of 67B, 67C, and 69G exerted significant effects on the motion modes of BRD9 and TAF1 (2). The BFEs of 67B, 67C, and 69G to BRD9 and TAF1(2) predicated by the MM-GBSA approach indicated that alterations in the binding enthalpy due to the inhibitors' associations with BRD9 compared with those with TAF1(2) were are primarily responsible for the binding selectivity of ligands to BRD9 over TAF1 (2). The results revealed that the enthalpy changes played critical roles in the selectivity recognition of ligands toward BRD9 and TAF1(2), which indicated that 67B and 67C could more favorably bind to TAF1(2) than BRD9, while 69G had better selectivity toward BRD9 versus TAF1(2). The results obtained from the energy contributions of individual residues indicated that three common residues, namely (I53 and F1536), (A54 and V1537), and (A96 and S1579) in (BRD9 and TAF1 (2)), generated significant differences in the associations of 67B, 67C, and 69G with BRD9 and TAF1 (2), indicating that these residues could be considered as hot spots for designing effectively selective inhibitors toward BRD9 and TAF1 (2). We also expect that this work could aid in obtaining a deeper understanding of the selectivity mechanism of inhibitors and provide theoretical guidance for the design of novel selective inhibitors targeting BRD9 and TAF1(2).